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In this paper recent substantial progress in applying the density-matrix renormalization- 
group (DMRG) to the simulation of the time-evolution of strongly correlated quantum systems 
in one dimension is reviewed. Various approaches to generating a time-evolution numerically are 
considered. The key focus of this review is on current strategies to circumvent the limitations 
of the quite small subspace well approximated by DMRG, by either enlarging or changing it 
as time evolves. All these approaches can be extended to the simulation of mixed, i.e. finite- 
temperature states. While this paper is phrased almost entirely in standard DMRG language, 
I finish by considering the alternative formulation of time-evolutions in the language of matrix 
product states which is less well-known but conceptually more powerful. 
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1. Introduction 

More than ten years after the invention of the 
Density-Matrix Renormalization Group (DMRG) by 
Steve White, 1,2 this method has become the method of 
choice for the numerical simulation of the equilibrium 
properties of strongly correlated one-dimensional quan- 
tum systems. 3,4 Ground-state properties of both bosonic 
and fermionic systems have been calculated often at al- 
most machine precision and comparatively low compu- 
tational cost. While these results have been of prime im- 
portance in understanding the details of gapped Haldane 
systems or critical Luttinger liquids, to name but a few 
applications, for a long time hardly any attention had 
been paid to the time-evolution of strongly correlated 
quantum systems, both due to the comparative lack of 
experimental input in the past and to the inherent diffi- 
culties of actually calculating such time-evolutions. The 
last years have seen an increasing number of experimen- 
tal results on non-trivial time-evolutions. Perhaps most 
spectacular was the recent mastery of storing ultracold 
bosonic atoms in a magnetic trap superimposed by an 
optical lattice: This has allowed to drive these atoms, at 
will, by time-dependent variations of the optical lattice 
strength, from the superfluid (metallic) to the Mott insu- 
lating regime. These regimes are linked by one of the key 
phase transitions in strongly correlated systems. 5 Quite 
generally, progress in the fields of nanoelectronics and 
spintronics raises the question how (strongly correlated) 
quantum many-body systems react to external time- 
dependent perturbations and how transport can be cal- 
culated quantitatively also far from the linear-response 
regime. On the computational physicists' side, the last 
twelve months or so have seen an extraordinary surge of 
activity in developing DMRG variants applicable to time- 
evolutions, such that by now we have various very pow- 
erful and conceptually innovative DMRG algorithms for 
pure as well as mixed quantum states, for non-dissipative 
as well as dissipative dynamics. 

The relationship between these various algorithms is 



often extremely close, and it is not always immediately 
obvious where the differences are of a fundamental or of 
a superficial nature. My goal in this paper is to present 
these algorithms ahistorically and to proceed as follows: 
after outlining the fundamental questions and difficulties, 
I first show that all mixed states (i.e. density matrices) 
can be cast into the form of pure states such that the 
entire discussion of this review can focus on the time- 
evolution of pure states. I move on to discussing various 
ways of actually generating time-evolution generated by 
Hamiltonians living in state spaces by far too large for 
exact diagonalization. This discussion is independent of 
DMRG specifics and serves to illustrate that in many of 
the proposals there is some freedom as to the choice of 
the actual time-evolution generator. DMRG takes center 
stage in the next section where various DMRG strate- 
gies of finding suitable state (sub)spaces to operate on 
are discussed: for all systems of interest here, the full 
state space is by far too large for numerical treatment. 
Recently, it has been shown that essentially all of the 
above can be expressed in a conceptually very beautiful 
form in the language of matrix product states and that 
computational performance and precision might be in- 
creased that way. I will attempt to clarify the connection 
between matrix product state and DMRG time-evolution 
algorithms. Throughout this paper, I assume familiar- 
ity with the basics of infinite-system and finite-system 
DMRG. 2-4 For notational brevity, I will consider time- 
independent Hamiltonians in the formulae; almost all re- 
sults carry over trivially to time- dependent Hamiltonians 
except where mentioned. 

Essentially all physical quantities of interest involving 
time can be reduced to the calculation of either equal- 
time n-point correlators such as the (1-point) density 

(«i(t)) = (m\ni\m) = (v^w-^'iv) (i) 

or unequal-time n-point correlators such as the (2-point) 
real-time Green's function 

G tJ (t) = M4(t)c,(0)M = Me+^cte-^cM). (2) 
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This expression can be cast in a form very close to Eq. (1) 
by introducing \4>) — Cj\tp) such that the desired correla- 
tor is then simply given as an equal-time matrix element 
between two time-evolved states, 



G ij (t) = m)\4m))- 



(3) 



If both \ip(t)) and \4>(t)} can be calculated, a very appeal- 
ing feature of this approach is that Gij(t) can be eval- 
uated in a single calculation for all i and t as time pro- 
ceeds. Frequency-momentum space is then reached by a 
double Fourier transformation. Obviously, finite system- 
sizes and edge effects as well as algorithmic constraints 
will impose physical constraints on the largest times and 
distances \i — j\ or minimal frequency and wave vectors 
resolutions accessible. Nevertheless, this approach might 
emerge as a very attractive alternative to the current 
very time-consuming calculations of G(k,uj) using the 
dynamical DMRG. 6 - 7 

The fundamental difficulty of obtaining the above cor- 
relators become obvious if we examine the time-evolution 
of the quantum state \ij){t — 0)) under the action of 
some time-independent Hamiltonian H \ip n ) = E n \ip n ). If 
the eigenstates \ip n ) are known, expanding \ip(t = 0)) = 
J2 n c n \ip n ) leads to the well-known time evolution 



(4) 



where the modulus of the expansion coefficients of \ip(t)) 
is time-independent. Now normally we are forced to 
restrict ourselves to some smaller state space, as rel- 
evant Hilbert space dimensions exceed computational 
resources. A sensible Hilbert space truncation is given 
by a projection onto the large-modulus eigenstates. In 
strongly correlated systems, however, we usually have 
no good knowledge of the eigenstates. Instead, one uses 
some orthonormal basis with unknown eigenbasis expan- 
sion, \m) — J2 n Q-mn\ipn}- The time evolution of the state 
= °)> = E m d ™(0)M thcn rcads 

M*)> = E (H d ™(0K m e- iiJ "*j |m) = J2d m (t)\m), 

m \ n / m 

(5) 

where the modulus of the expansion coefficients d m (t) is 
time- dependent. For a general orthonormal basis, Hilbcrt 
space truncation at one fixed time (i.e. t = 0) will there- 
fore not ensure a reliable approximation of the time evo- 
lution. Also, energy differences matter in time evolution 
due to the phase factors e~' 1 ^ En ~ E ^ t in \d m (t)\ 2 . Thus, 
the sometimes justified hope that DMRG yields a good 
approximation to the low-energy Hamiltonian is of lim- 
ited use. 

2. Time-evolution of mixed states 

After the previous discussion on the difficulties of sim- 
ulating the time-evolution of pure states in subsets of 
large Hilbert spaces it may seem that the time-evolution 
of mixed states (density matrices) is completely out of 
reach. It is however easy to see that a thermal density 
matrix pp = exp[— /3H) can be constructed as a pure 
state in an enlarged Hilbert space and that Hamiltonian 



dynamics of the density matrix can be calculated con- 
sidering just this pure state (dissipative dynamics being 
more complicated). In the DMRG context, this has first 
been pointed out by Verstraete, Garcia-Ripoll and Cirac 8 
and Zwolak and Vidal, 9 using essentially information- 
theoretical language; it has also been used previously in 
pure statistical physics language in e.g. high-temperature 
series expansions. 10 

To this end, consider the completely mixed state p = 
I. Let us assume that the dimension of the local physical 
state space {|o"i)} of a physical site is N. Introduce now a 
local auxiliary state space {In}} of the same dimension 
N on an auxiliary site. The local physical site is thus 
replaced by a rung of two sites, and a one-dimensional 
chain by a two-leg ladder of physical and auxiliary sites 
on top and bottom rungs. Prepare now each rung i in 
the Bell state 



N 



E Wi T i) 



(6) 



Other choices of |f/>o) are equally feasible, as long as they 
maintain in their product states maximal entanglement 
between physical states \afj and auxiliary states |ri). 
Evaluating now the expectation value of some local op- 
erator O l a acting on the physical state space with respect 
to IV'o), one finds 

moi\^) = E E ' 



-n <j'=t' 



TV 



The double sum collapses to 



moi\%) 



N 



N ^ 



and we see that the expectation value of % a with respect 
to the pure state living on the product of physical 
and auxiliary space is identical to the expectation value 
of % a with respect to the completely mixed local physical 
state, or 



where 



(01) = Tr^Ol 



pi = Tr T \r )(%\- 



(7) 



(8) 



This generalizes from rung to ladder using the density 
operator 



where 



/So = Tt^lVoX^ol 



m = n i^o> 

i=l 



(9) 



(10) 



is the product of all local Bell states, and the conver- 
sion from ficticious pure state to physical mixed state is 
achieved by tracing out all auxiliary degrees of freedom. 
At finite temperatures [3 > one uses 



P(3 



e -f3H/2 Ie - H/2 = Tr r e-^/ 2 |V>o><Vo|e- 



0H/2 
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where we have used Eq. (9) and the observation that the 
trace can be pulled out as it acts on the auxiliary space 
and e -P H / 2 on the real space. Hence, 



p fj = Tr r |V>/3)(#|, 



(11) 



where = e l3H / 2 \tp )- Similarly, this finite- 

temperature density matrix can now be evolved in 
time by considering \tpp(t)) = e~ lHt \tpp(0)) and 
Pp(t) = Tr T \il)p(t)){ipp(t)\. The calculation of the finite- 
temperature time-dependent properties of, say, a Hub- 
bard chain, therefore corresponds to the imaginary-time 
and real-time evolution of a Hubbard ladder prepared to 
be in a product of special rung states. Time evolutions 
generated by Hamiltonians act on the physical leg of the 
ladder only. As for the evaluation of expectation values 
both local and auxiliary degrees of freedom are traced 
on the same footing, the distinction can be completely 
dropped but for the time-evolution itself. Code-reusage 
is thus almost trivial. Note also that the initial infinite- 
temperature pure state needs only M = 1 block states 
to be described exactly in DMRG as it is a product state 
of single local states. Imaginary-time evolution (lowering 
the temperature) will introduce entanglement such that 
to maintain some desired DMRG precision M will have 
to be increased. 

3. Generating time-evolutions in large state 
spaces 

To generate time-evolutions in large state spaces, sev- 
eral approaches are feasible to integrate the Schrodingcr 
equation 



4 t \m) 



(12) 



Here, E = (ip(0)\H\ip(0)} . If H itself is time-independent, 
energy is conserved and for reasons of numerical effi- 
ciency a trivial phase factor of exp(— iEt) may be ex- 
tracted from \ip(t)). In the following I assume energy 
shifts to set E = 0. 

On the one hand, one may now directly integrate the 
differential equation forward in time, choosing some in- 
finitesimal step size At. To lowest order in At, the solu- 
tion is 



m+At)) = (i-iH(t)M)m))- 



(13) 



Many improvements on this simple-minded approach 
exist. In the fourth-order Runge-Kutta algorithm, 11-14 
which is perhaps the most popular, one calculates, start- 
ing from \ip(t)), the following four Runge-Kutta vectors, 



\ki) = 


-iAtHm)), 


(14) 


1*2) = 


-iAtH{\m) + \\ki)), 


(15) 


1*3) = 


-\AtH{m)) + ^2», 


(16) 


\k 4 ) = 


-\AtH(m)) + \h)). 


(17) 



The wave function \ip(t + At)} is then given by 

W + At)) = |V>(i)} + ^|fcl} + ^|fc2> + ^3> + ^|fc4}. (18) 



Note that numerical efficiency hinges centrally on the 
fast evaluation of matrix- vector products H\t(j(t)). 

It is important to note that Runge-Kutta does not pre- 
serve unitarity. This can be improved by using e.g. the 
unitary Crank-Nicholson time-evolution, 13 in its most 
simple version given as 



m+At)) 



1 - iH{t)At/2 



(19) 



1 + iH{t)At/T 

The drawback of this approach is that the evaluation of 
the non-hermitian denominator is non-trivial. This may 
be done by considering the large linear equation system 

{l + \H{t)At/2)\4>) = m)), (20) 

such that 

\iP{t + At)) = {l-\H{t)At/2)\(f>). (21) 

The non-hermitian equation system (20) may be solved 
by using a biconjugate gradient approach. In fact, like 
all conjugate gradient approaches, its convergence can 
be greatly improved by preconditioning it, which here 
means the provision of the solution of some closely re- 
lated, exactly solvable equation system. Such a solution 
is trivially given here by considering the diagonal "equa- 
tion system" \<j>) = \ip(t)), as At — ► ensures that Eq. 
(20) is diagonally dominated. 

Instead of considering time-evolution as a differen- 
tial equation to be integrated numerically, one may also 
consider 15 for time-independent Hamiltonians the for- 
mal solution of Eq. (12), the time-evolution operator 
exp(— iHt). This avoids the numerical delicacies of the 
integration of the Schrodinger equation, but introduces 
new ones. It is known that the numerical evaluation of 
matrix exponentials is non-trivial; 16 for example, it is 
not a good strategy to consider their Taylor expansion. 
It seems that among the most efficient approaches is 
the so-called Krylov subspace approximation. Remem- 
bering that we are never interested in exp(— iHt), but 
in exp(— iHt)\tjj), one may exploit that in DMRG H\ip) 
must be available efficiently, and form the Krylov space 
by successive Schmidt-Gram orthonormalization of the 
set {\tp}, —iHt\ip), (-iHt) 2 \ip), . . .} assumed normal- 
ized). As in the closely related Lanczos algorithm, the 
n Krylov vectors thus obtained form the orthonormal 
columns of a matrix V such that —iHt is approximated 
regarding its extreme eigenvalues by VTV T , where T is 
an n x n tridiagonal matrix, to excellent precision even 
for relatively small n. The exponential is then given by 
the first column of V^expT, where the latter exponential 
is now much easier to calculate. 

Another way of evaluating the above matrix exponen- 
tial that has its origin in Quantum Monte Carlo meth- 
ods is based on the Suzuki- Trotter decomposition. 17 ' 18 
Its main usefulness is for Hamiltonians with nearest- 
neighbor interactions. In that case we can separate inter- 
actions on odd bonds from interactions on even bonds, 
such that the decomposition reads H = H\ + H 2 ; Hi = 



N/2 



Here, hi is the local 



Hamiltonian linking sites i and i + 1. As neighbouring 
local Hamiltonians will in general not commute, we have 
[Hi,H 2 ] ^ 0. However, all terms in Hi or H 2 commute. 
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The first-order Trotter decomposition of the infinitesimal 
time-evolution operator then reads 

exp(-LHAi) = exp(-LHiAf) exp(-i£T 2 Ai) + 0(At 2 ). 

(22) 

The second-order Trotter decomposition reads 

exp(-LHAt) = (23) 

p -iH 1 At/2 e -iH 2 At e -iH 1 At/2 + Q(At 3 ). 

A frequently used forth-order Trotter decomposition 
reads 



exp(-LHAt) = 



(24) 



n 



e -ipiH! At/2 e -i P iH 2 At e -ipiH 1 At/2 



+ 0{At 5 ), 



where all pi = 1/(4 — 4 1 / 3 ), except p 3 = 1 — 4pi < 0, 
corresponding to evolution backward in time. 

For a nth-order Trotter decomposition, the error made 
in one time step At is of order At n+1 . To reach a given 
time t one has to perform t/ At time-steps, such that the 
error grows (at worst) linearly in time t and the resulting 
error is bounded by an expression of order (At) n t. 

4. DMRG state space selection under time- 
evolution 

All time-evolution schemes for DMRG so far can be 
classified into three categories as regards their state space 
selection. Static Hilbert space DMRG methods keep the 
truncated Hilbert space found optimal for \i[)(t = 0)). 
What I will call dynamic Hilbert space DMRG methods 
try to enlarge the truncated Hilbert space starting from 
the one optimal to approximate \ip(t — 0)), such that it is 
big enough to accommodate (i.e. maintain large overlap 
with the exact result) \tp{t)} for a sufficiently long time 
to a very good approximation. More recently, adaptive 
Hilbert space DMRG methods keep the size of the trun- 
cated Hilbert space fixed, but try to change it as time 
evolves such that it also accommodates \ip(t)) to a very 
good approximation. 

4-1 Static time- dependent DMRG 

Cazalilla and Marston 11 were the first to exploit 
DMRG to systematically calculate time-dependent quan- 
tum many-body effects. They studied a time-dependent 
Hamiltonian H(t) = H(0) + V{t), where V(t) encodes 
the time-dependent part of the Hamiltonian. After ap- 
plying a standard DMRG calculation to the Hamilto- 
nian H(t = 0), the time-dependent Schrodinger equation 
was numerically integrated forward in time. The effec- 
tive Hamiltonian in the reduced Hilbert space was built 
as H e e(t) = H e s(0) + V e s(t), where i? e ff(0) was taken 
as the last superblock Hamiltonian approximating H(0). 
V e g(t) as an approximation to V was built using the rep- 
resentations of operators in the final block bases. The ini- 
tial condition was obviously to take IV'(O)) as the ground 
state obtained by the preliminary DMRG run. This pro- 
cedure amounts to fixing the reduced Hilbert space at 
that optimal at t = 0, and projecting all wave functions 
and operators onto it. 



As an application, Cazalilla and Marston have consid- 
ered a quantum dot weakly coupled to noninteracting 
leads of spinless fermions, where time-dependency is in- 
troduced through a time-dependent chemical potential 
coupling to the number of particles left and right of the 
dot, 



V(t) = -5n H (t)N R - Sfi L (t)N L , 



(25) 



which is switched on smoothly at t = 0. Setting 5/j.l — 
—SfiR, one may gauge-transform this chemical potential 
away into a time-dependent complex hopping from and 
to the dot, 



t q (t) = t q exp 



f 

J — < 



dt'5fj, L (t') 



(26) 



The current is then given by evaluating the imaginary 
part of the local hopping expectation value. Obviously, in 
a finite system currents will not stay at some steady state 
value, but go to zero on a time scale of the inverse system 
size, when lead depletion has occurred. The calculation 
was helped by the fact that the time-dependent coupling 
is in the center of the system which is modelled exactly 
by explicit sites in DMRG. 

In this approach the hope is that an effective Hamilto- 
nian obtained by targeting the ground state of the t = 
Hamiltonian is capable to catch the states that will be 
visited by the time-dependent Hamiltonian during time 
evolution. This approach must however break down after 
relatively short times as the full Hilbert space is explored. 

4-2 Dynamic time- dependent DMRG 

Several attempts have been made to enlarge the re- 
duced Hilbert space using information on the time- 
evolution, such that the time-evolving state has large 
support on that Hilbert space for longer times. Whatever 
procedure for enlargement is used, the problem remains 
that the number of DMRG states M grows with the de- 
sired simulation time as they have to encode more and 
more different physical states. As calculation time scales 
as M 3 , this type of approach will meet its limitations 
somewhat later in time. 

All enlargement procedures rest on the ability of 
DMRG to describe - at some numerical expense - small 
sets of states ("target states") very well instead of just 
one. This is achieved by considering for the DMRG state 
selection reduced density matrices formed as sums of 
some suitably chosen pure state projectors, 



Tr, 



p = Tr £ ^a i |^)(^|, (27) 



= 1 for normalized Tr# is the standard 

DMRG trace over the environment. 

The simplest approach 13 is to target the set 
{|Vi>} = {|V(0)),i?|V(0)),i? 2 |^(0)),...}. Alternatively, 
one might consider the Krylov vectors formed from this 
set. Results improve, but not decisively. 

A much more time-consuming, but also much more 
performing approach has been demonstrated by Luo, Xi- 
ang and Wang. 12 They use a density matrix that is given 
by a superposition of states \ip(U)) at various times of 
the evolution, p = Sil'o a *IV'(^))(V'(*i)l wi th a i = 1 
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for the determination of the reduced Hilbert space. Of 
course, these states are not known initially; it was pro- 
posed by them to start within the framework of infinite- 
system DMRG from a small DMRG system and evolve it 
in time. For a very small system this procedure is exact. 
For this system size, the state vectors \i/)(U)) are used 
to form the density matrix. This density matrix then de- 
termines the reduced Hilbert space for the next larger 
system, taking into account how time-evolution explores 
the Hilbert space for the smaller system. One then moves 
on to the next larger DMRG system where the procedure 
is repeated. This is of course very time-consuming, and 
it would also be of advantage to extend the procedure to 
the usage of the finite-system DMRG method. 

Schmitteckert 15 has computed the transport through 
a small interacting nanostructure using an Hilbert space 
enlarging approach, based on the time evolution opera- 
tor. To this end, he splits the problem into two parts: 
By obtaining a relatively large number of low-lying 
eigenstates exactly (within time-independent DMRG 
precision), one can calculate their time evolution ex- 
actly. For the subspace orthogonal to these eigenstates, 
he implements the matrix exponential \tp(t + At)) = 
exp(—iHAt)\tp(t)) using the Krylov subspace approxi- 
mation. For any block-site configuration during sweep- 
ing, he evolves the state in time, obtaining \ip(U)) at 
fixed times t%. These are targeted in the density matrix, 
such that upon sweeping forth and back a Hilbert space 
suitable to describe all of them at good precision should 
be obtained. For numerical efficiency he carries out this 
procedure to convergence for some small time, which is 
then increased upon sweeping, bringing more and more 
states \tp(ti)) into the density matrix. 

4-3 Adaptive time- dependent DMRG 

All approaches mentioned so far imply that to main- 
tain good precision over time the size of the reduced 
Hilbert space has to grow with desired simulation time, 
imposing huge algorithmic cost. Time-dependent DMRG 
using adaptive Hilbert spaces, that "follow" the state to 
be approximated optimally, has first been proposed in 
essentially identical form by Daley Kollath, Schollwock 
and Vidal 13 and White and Feiguin. 14 Both approaches 
are efficient implementations of an algorithm for the clas- 
sical simulation of the time evolution of weakly entangled 
quantum states invented by Vidal 19 ' 20 [time-evolving 
block-decimation (TEBD) algorithm]. The TEBD algo- 
rithm was originally formulated in the matrix product 
state language. For simplicity, I will explain the algo- 
rithm in its DMRG context; a detailed discussion of the 
very strong connection between adaptive time-dependent 
DMRG and the original simulation algorithm is given by 
Daley et al.: 13 it turns out that DMRG naturally at- 
taches good quantum numbers to state spaces used by 
the TEBD algorithm, allowing for the usual drastic in- 
creases in performance due to the use of good quantum 
numbers. 

For each block-site configuration (system-sites- 
environment: S»»E) during DMRG calculations, the 
truncation is always carried out at the position of the 
sites; hence it can react to changes of the physical state 



there. On the other hand, as these sites are explicit, 
we can carry out time-evolution on this bond even 
exactly. Time evolution in the adaptive time-dependent 
DMRG is therefore generated using a Suzuki- Trotter 
decomposition: 18 Expanding Hi and H 2 into the local 
Hamiltonians, one infinitesimal time step t — > t + At 
may be carried out (taking the most simple case of a 
first-order Suzuki- Trotter decomposition) by performing 
the local time-evolution on (say) all even bonds first 
and then all odd bonds. We may therefore carry out 
an exact time-evolution by performing one finite-system 
sweep forward and backward through an entire one- 
dimensional chain, with time-evolutions on all even 
bonds on the forward sweep and all odd bonds on the 
backward sweep, at the price of the Trotter error. This 
procedure necessitates that \ip(t)) is available in the 
right block bases for the current block-site configura- 
tion, which is ensured by carrying out the reduced basis 
transformations on \ij){t)) that in standard DMRG form 
the basis of White's prediction method. 21 The decisive 
idea of Vidal 19,20 in the TEBD algorithm was now to 
carry out a new Schmidt decomposition (in DMRG 
language: reduced density matrix formation) and make 
a new choice of the most relevant block basis states for 
\4>(t)) after each local bond update. Therefore, as the 
quantum state changes in time, so do the block bases 
such that an optimal representation of the time-evolving 
state is ensured. Choosing new block bases changes 
the effective Hilbert space, hence the name adaptive 
time-dependent DMRG. 

An appealing feature of this algorithm is that it can be 
very easily implemented in existing finite-system DMRG. 
One uses standard finite-system DMRG to generate a 
high-precision initial state \ip(0)) and continues to run 
finite-system sweeps, one for each infinitesimal time step, 
merely replacing the large sparse-matrix diagonalization 
at each step of the sweep by local bond updates for the 
odd and even bonds, respectively. 

An extensive error discussion has been carried out by 
Gobert et al. 23 Two main sources of errors occur, the 
Trotter error due to the Suzuki- Trotter decomposition, 
and the truncation error due to the loss of information 
about the quantum state at each DMRG projection onto 
a reduced basis. As discussed above, the Trotter error 
scales with time as 0(At n t) for a nth order Suzuki- 
Trotter decomposition. Of course, there are prefactors 
associated depending on the details of the decomposi- 
tion. In the setup of Gobert et al. the error scales (max- 
imally) linearly with system size L, and overall it is thus 
of order (At) n Lt. Further errors are produced by the 
truncations (reduced basis transformations) after each 
bond update during an infinitesimal time step. While 
the truncation error e that sets the scale of the error of 
the wave function and operators is typically very small 
in DMRG, here it will strongly accumulate over time as 
O(LtfAt) truncations are carried out up to time t. This 
is because the truncated DMRG wave function has norm 
less than one and is renormalized at each truncation by 
a factor of (1 — e) _1 > f . Truncation errors should there- 
fore accumulate roughly exponentially with an exponent 
of eLt/At, such that eventually the adaptive t-DMRG 
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will break down at too long times. In practice, satura- 
tion effects in errors and partial compensations of errors 
in observables will lead to somewhat slower error growth. 
The very different growth behaviour of these errors leads 
to a well-defined "runaway time" , at which the Trotter 
error is replaced by the accumulated truncation error as 
dominant error. It can be seen in a sharp onset of error 
growth for physical observables, either be comparison to 
higher-precision (larger M) calculations or to exact solu- 
tions, when available. The errors mentioned can be well 
controlled by increasing M, the runaway time, which sig- 
nals that adaptive t-DMRG will break down in the im- 
minent future, growing linearly 23 with M. 

Applications of adaptive time-dependent DMRG so 
far for the (time-dependent) Bose-Hubbard model, 13 ' 22 
spin-1 Heisenberg chains, 14 and far-from-equilibrium 
states in spin-1/2 Heisenberg chains 23 have demon- 
strated that it allows to access large time scales (several 
hundred in units of inverse interaction energy or band- 
width) very reliably. 

Adaptive t-DMRG based on Suzuki- Trotter decompo- 
sitions is for all practical purposes restricted to nearest- 
neighbour interactions. In a second adaptive approach 
without that limitation, White and Feiguin 24 ' 25 exploit 
that there is no conceptual need in DMRG to combine 
the actual time-evolution and the adaptation of the state 
spaces into one step. In fact, even the time-scales St on 
which time-evolution is discretized and At on which state 
spaces are adapted need not be identical: the empirical 
observation that even a static state space remains a good 
choice for some finite time, allows to choose At ^> St. 
The actual time-evolution may be carried out using any 
of the non- Trotter methods described above; White and 
Feiguin have chosen fourth-order Runge-Kutta integra- 
tion, which is a valid and convenient, but perhaps not the 
ultimate choice. One may alternatively use, for example, 
Crank- Nicholson integration or Krylov- vector based ex- 
ponentials. 26 ' 27 

To adapt the state space, at intervals At, an analy- 
sis is carried out on the states that will appear in the 
near future. To that purpose, several DMRG sweeps are 
carried out at a fixed time t. For each block-site config- 
uration during these sweeps, a Runge-Kutta integration 
up to t + At is carried out. The reduced density ma- 
trix that during sweeping controls the choice of the re- 
duced block bases is formed from the present state and 
the states representing the near future. Here resides an 
ambiguity. While the Runge-Kutta vectors encode 
time-evolution up to t + At, the density matrix choice 
p s = Tr E {\^{t))(^(t)\ + \h)(h\) is not found to be 
optimal. The Runge-Kutta vectors are essentially encod- 
ing derivatives of state vectors, while what one wants to 
approximate are state vectors at times up to t + At. The 
density matrix during sweeping is therefore constructed 
as 

4 

Ps = 1^X>(*0>W*i)l (28) 

i=l 

where t 1 =t,t 2 = t + At/3, t 3 = t + 2At/3, t 4 = t + At. 
Here, of course modifications are conceivable. The state 
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vectors at the intermediate times t + At/3 and t + A2t/3 
are given by 

m+At/3)) = m))+ (29) 

^ [31|fc!> + U(\k 2 ) + |fc 3 » - 5|fc 4 >] 

|V(i + A2t/3)) = |V(*)> + (30) 

^-[16|fc 1 )+20(|fc 2 ) + |fc 3 ))-2|fc 4 )] 

They can be calculated from exploiting that Runge- 
Kutta methods are exact for polynomials up to some de- 
gree, which establishes a connection between the Runge- 
Kutta vectors and the state at all intermediate times. 

This approach is free of the Trotter error, but the error 
analysis of the truncation error should remain valid. In 
fact, there might be for long-time simulations additional 
errors from finding the approximate Hamiltonian (i.e. the 
reduced state space on which the full Hamiltonian is pro- 
jected) by carrying out a time-evolution using just that 
approximate Hamiltonian. That there is some such error 
can be seen from the observation that for At — > oo one 
essentially recovers static time-dependent DMRG with 
its limitations. Of course, for any finite At, the adap- 
tive DMRG will be (much) better than static DMRG. 
At therefore has to be chosen optimally: small enough to 
update the Hilbert space often enough, large enough to 
avoid too much numerical overhead. At the time of writ- 
ing, a complete error analysis still remains to be done. 
The price to be paid for the elimination of the Trotter er- 
ror and the generalization to longer-ranged interactions 
is that during the finite-system sweeps at intervals of At, 
a substantial number of time-consuming H\ip) multipli- 
cations have to be carried out, which can be completely 
avoided in the Trotter-based approach. 

As the time-evolution of density matrices can be simu- 
lated by using pure states in enlarged state spaces, all the 
above carries over to the simulation of imaginary time- 
evolutions. Two comments are in order, (i) The com- 
pletely mixed state can be described in DMRG using 
M = 1. After a very long imaginary time-evolution, it 
will reach the ground state, for which - given years of 
ground state DMRG experience - substantial values of 
M (typically hundreds) have to be chosen for high pre- 
cision. It is therefore a sensible approach to carry out 
imaginary time evolutions fixing some truncation error 
and to choose M at each time step such that it is not 
exceeded. This will amount to keeping more and more 
states as temperature is being lowered, (ii) Imaginary- 
time evolution is not unitary, but reduces norms. This 
leads to compensating effects for the accumulation of er- 
rors and imaginary-time evolution is numerically easier 
to handle than real-time evolution. 

5. An alternative formulation: matrix-product 
states 

As has been pointed out by various authors, 28-32 the 
ansatz states that DMRG generates are very closely re- 
lated to position- dependent matrix-product states (also 
known as finitely correlated states). 33 ' 34 However, there 
are subtle, but crucial differences between DMRG states 
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and matrix-product states. 32 ' 35 Recently, consistent use 
of matrix-product states has given rise to new algorithms 
closely related to DMRG 8,35 that are of high relevance 
for time-evolution. Let me first state some properties of 
matrix-product states. Matrix-product states are simple 
generalizations of product states of local states, which 
we take to be on a chain, 



|er) = ® |<7 2 ) <X> . . . <X> \<Jl), 



(31) 



obtained by introducing linear operators .Aj[<7i] depend- 
ing on the local state. These operators map from some 
M-dimensional auxiliary state space spanned by an or- 
thonormal basis {|/3)} to another M-dimensional auxil- 
iary state space spanned by {|a)}: 



M [a i 



a/3 



(AiN) a/3 |a)(/3|. 



(32) 



One may visualize the auxiliary state spaces to be located 
on the bonds + 1) and (i — The operators are 
thus represented by M x M matrices (Ai[ai]) a p; M will 
be seen later to be the number of block states in DMRG. 
We further demand for reasons explained below that 



]A i [a i ]Al[a i ]=I. 



(33) 



A position-dependent unnormalized matrix-product 
state for a one-dimensional system of size L is then given 

by 



IV>) = E [(HIIIMvMr) k>, 



(34) 



where (</>l| and \<pn) are left and right boundary states 
in the auxiliary state spaces located in the above visual- 
ization to the left of the first and to the right of the last 
site. They are used to obtain scalar coefficients. Position- 
independent matrix-product states are obtained by mak- 
ing Eq. (32) position-independent, Aj[<Tj] — > A[<Ti\. For 
simplicity, we shall consider only those in the following. 

For periodic boundary conditions, boundary states are 
replaced by tracing the matrix-product: 



E Tr 

{er} 



.i=l 



\a). 



(35) 



The best-known matrix-product state is the valcncc- 
bond-solid ground state of the bilinear-biquadratic S — 
1 Afflcck-Kcnnedy-Licb-Tasaki Hamiltonian, 36,37 where 
M = 2. 

If we consider operators Oj and Oj + i acting on sites 
j and j + I, for a periodic boundary condition matrix- 
product state the correlator C(l) — (ip\OjOj + i\if>) / (ip\ip) 
is given by 38 



C(l) 



TvO/~ 



TrI 



(36) 



where we have used the following mapping 30,38 from an 
operator O acting on the local state space to a M 2 - 
dimensional operator O acting on products of auxiliary 
states |/3/3') = \/3) <g> \/3'): 



O 



E 



{a'\6\a)A*[a'\® A[a\. 



(37) 



Note that A* stands for A complex-conjugated only as 
opposed to A'. Evaluating Eq. (36) in the eigenbasis of 
the mapped identity, I, we find that in the thermody- 
namic limit L 



oo 



M / \ \ l 



(38) 



with £j = — 1/ In |Aj|. The Aj are the eigenvalues of I, and 
the Cj depend on O. This expression holds because due to 
Eq. (33) |Aj| < 1 and Ai = 1 for the eigenstate (/3/3'|Ai> = 
S/3/3'. Equation (33) is thus seen to ensure normalizability 
of matrix product states in the thermodynamic limit. 
Generally, all correlations in matrix-product states are 
either long-ranged or purely exponentially decaying. 

That a DMRG calculation retaining M block states 
produces states closely related toMxM matrix-product 
states as introduced above was first realized by Ostlund 
and Rommer. 28 Let us consider (in the framework of 
finite-system DMRG) the reduced basis transformation 
to obtain the block states |mj) of a right block including 
sites i,. . . ,L from the right block states |mj+i) of sites 
i + 1, . . . , L and the states |<7j) of site i: 

(m i+1 ai\mi) = (Ai) mi . mi+irTi = (i,[(r,]) mi;m , +1 , (39) 

such that 

\rm)= E ( A i[(Ti]) mi ;m i+1 \(Ti) O \m i+1 ) . (40) 

The reduced basis transformation matrices ^[ct,] auto- 
matically obey Eq. (33), which here ensures that {|mi)} 
is an orthonormal basis provided {|mj + i)} is one, too. 
For the left block of sites 1, . . . , i, similarly define 

(mi-i<ri\mi) = (Ai) mi _ iai -„ H = (i,[ff,]) mi _ i;m „ (41) 

which implies 

Y,MWi]M<n]=*, (42) 

a variation of Eq. (33) that also ensures normalization. 

We may now use Eq. (40) for a backward recursion to 
express |mj+i) via |mj+2) and so forth. I will show re- 
sults for the right blocks, the left blocks can be treated 
analogously. There is a (conceptually irrelevant) compli- 
cation as the number of block states for very short blocks 
is less than M. For simplicity, I assume that N N = M, 
and stop the recursion at the shortest right block of size 
N that has M states, such that 



i m *> = E E 

er BB oi a. 



(Ai[ai] . . . A L _f f [a L _f J ]) mi (T RB x 



\<7i...a L _fi}®\<rR B ), (43) 

where we have boundary-site states \(Trb) = 
\ a L-N+i hence 

K> = E (M (J t]---A L _f f [(7 L ^^}) m ^( TRB \^...(7L}- 

<7i,...,a L 

(44) 

DMRG thus generates position-dependent M x M 
matrix-product states as block states for a reduced 
Hilbert space of M states; the auxiliary state space to 
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a local state space is given by the Hilbert space of the 
block to which the local site is the latest attachment. 
DMRG looks for the ground state of the full chain as 
the variational optimum in energy in a space spanned 
by products of two local states and two matrix-product 
states, 



2 \mi-i<Jiai + im,i + 2) 



(45) 



Encoding the ^-coefficients as a set of M x M ma- 
trices parametrized by a^i+i, * TOi _ imi+2 [<Ji<T i+1 ] = 
VWiow+im,^: we have 

IV>> = X](^iv+>iv+i]---^-i[^-i]*[<W+i] x 
{"■} 

A i+2 [a i+ 2] . . . A L _f l [a L _ 1 j-})cr LB ,(T RB 
\ai...a L ). (46) 

The effect of the finite-system DMRG algorithm can 
be seen from Eq. (46) to be a sequence of local optimiza- 
tion steps of the wave function that have two effects: on 
the one hand, the variational coefficients vf^crjCj+i] are 
optimized, on the other hand, a new improved ansatz 
matrix A is obtained for the site added into the growing 
block, using the improved variational coefficients from ^> 
for the reduced density matrix to obtain A. 

On closer inspection, the DMRG state to describe \ip) 
is different from a true matrix-product state: There is an 
anomaly in that formal "translational invariance" of this 
state is broken by the indexing of by two sites. This 
suggests to consider instead a setup where * is indexed 
by one site: 

W) = E( i w+i^+i]-- J '-i[ ff -il*N x 

{<T} 

A i+1 [a i+1 ] . . ■ A L _ft[a L _ft])(r LB ,(T RB 
|<7!... o- L ). (47) 

In the block-site language of DMRG this would corre- 
spond to a S»E instead of an S»»E setup where bullets 
stand for sites. The (finite-system) DMRG algorithm for 
this S»E setup can now be reformulated in this picture 
as follows: sweeping forward and backward through the 
chain, one keeps for site i all A at other sites fixed and 
seeks the "J/j that minimizes the total energy. From this, 
one determines A, and moves to the next site, seeking 
ipi+i (or ipi-i, depending on the direction of the sweep), 
and so on, until all matrices have converged. 

The key point is that in the ansatz of Eq. (47) the 
variational optimum possible for a matrix product state 
of dimension M, whereas this is not the case for the 
ansatz of Eq. (46), albeit it may come very close. This is 
the major difference between DMRG ansatz and matrix 
product states. To see this consider how A[<Ji] and the 
variationally optimal ^[ai] are related in the new setup: 
A is obtained from a reduced density matrix which is 
formed from tracing out the environment in the pure 
state projector In the new setup the environ- 

ment is a block of dimension M. Linear algebra tells 
us that the maximum number of non-zero eigenvalues 



of the MiV-dimcnsional density matrix thus obtained 
is then bounded by M. "Truncation" down to the M 
highest- weight eigenstates of the reduced density matrix 
is hence an exact basis transformation of the variation- 
ally obtained state of minimal energy. There is no loss 
of information. Shifting the active site i through the sys- 
tem means taking out some old A, look for the \E" with 
minimal energy at its place, and find a new A from it. 
As there is a ^ corresponding to that old A (i.e. leading 
to the same energy of the state; it can be constructed 
explicitly), the minimal energy can only stay constant 
or go down, such that this sequence of improvements of 
the ansatz matrices is truly variational in the space of 
the states generated by the maps A and reaches a min- 
imum of energy within that space (there is of course no 
guarantee to reach the global minimum). By comparison, 
the setup S»»E leads to a reduced basis transformation 
with generic loss of information as the dimension of the 
environment is MN. It is hence not strictly variational. 

Indeed, in practical applications it has been observed 
that for translationally invariant systems with periodic 
boundary conditions the position dependency of ob- 
servables is not eliminated entirely during-finite sys- 
tem sweeps in standard DMRG. Dukclsky et al. 31 and 
Takasaki, Hikihara and Nishino 32 were the first to ob- 
serve that switching, after convergence is reached, from 
the S»»E scheme for the finite-system DMRG algorithm 
to a S»E scheme for further sweeps eliminates this de- 
pendency. 

Let us now turn back to the original topic of this 
paper. To calculate time-evolutions in matrix-product 
states, Verstraete, Garcia-Ri'poll and Cirac 8 consider the 
T = matrix product state of Eq. (35). The Ai[ai] 
arc now interpreted as maps from the tensor prod- 
uct of the two auxiliary states (dimension M 2 ) next 
to a site to the physical state space of the site of di- 
mension N. This can be written as: (Ai[a]) a p: Ai = 
E„, Z) Qi ft(^[ cr ])a/3l cr i)( Q: iAI- The \a) and \(3) are states 
of the auxiliary state spaces ai and bi on the left and right 
bonds emerging from site i. 

The state of Eq. (35) and its interpretation can be 
generalized to finite-temperature matrix product density 
operators 



p= E Tr 

{<T}{<7'} 



\*){*>\. 



(48) 



Here, the M\oio'^ now are maps from the tensor product 
of four auxiliary state spaces (two left, two right of the 
physical site) of dimension M 4 to the iV 2 -dimensional 
local density-operator state space. The most general case 
allowed by quantum mechanics for such mappings is for 
M to be a completely positive map. The dimension of M 
seems to be prohibitive for making practical use of this 
variational ansatz (48), but it is known from quantum 
information theory that M can be decomposed into a 
sum of d tensor products of A-maps as 



Miio-io-'j] = Mum] ® A \ 

o«=l 



(49) 



J. Phys. Soc. Jpn. 



Full Paper 



Ulrich Schollwock 9 



In general, d < NM 2 , but it will be important to realise 
that for thermal states d — N only. At T = 0, one re- 
covers the standard matrix product state, i.e. d = 1. In 
order to actually simulate p, Verstraete et al. 8 consider 
the purification 

L 



\4>. 



UPDO 



{er}{-r} 



JJ A i [a l T i ] 



(50) 



such that p = Ttt\iPmpdo){' i Pmpdo\- Here, ancilla state 
spaces {|Tj)} of dimension d have been introduced. 

In this form, a completely mixed state is obtained from 
matrices j4j[ (J i T i] cx l-S auTi , where M may be 1 and nor- 
malization has been ignored. This shows d = N. This 
state is now subjected to infinitesimal evolutions in imag- 



inary time, e 



-HAt 



up to the imaginary time (3/2. As they 



act on a only, the dimension of the ancilla state spaces 
need not be increased. Of course, for T = the state may 
also be efficiently prepared using standard methods. The 
introduction of an ancilla state space of dimension TV is 
the complete equivalent of the motivation of the auxiliary 
state space introduced in Sec. 2. 

The imaginary-time evolution is carried out after a 
Trotter decomposition into infinitesimal time steps on 
bonds. The local bond evolution operator Ui_i+i is con- 
veniently decomposed into a sum of djj tensor products 
of on-site evolution operators, 



Ui, i+1 = J2 u i 

k=l 

djj is typically small, say 2 to 4. Applying now the 
time evolution at, say, all odd bonds exactly, the aux- 
iliary state spaces are enlarged from dimension M to 
Mdjj. One now has to find the optimal approximation 
\tp(t + At)) to \iji(t + At)) using auxiliary state spaces 
of dimension M only. Hence, the state spaces of di- 
mension Mdjj must be truncated optimally to minimize 
|| + At)) - \tp(t + At))\\. If one uses the matrices 
composing the state at t as initial guess and keeps all 
^4-matrices but one fixed, one obtains a linear equation 
system for this A; sweeping through all ^4-matrices sev- 
eral times is sufficient to reach the fixed point which is 
the variational optimum. As temperature is lowered, M 
will be increased to maintain a desired precision. Once 
the thermal state is constructed, real-time evolutions 
governed by Hamiltonians can be calculated similarly. 
It therefore turns out that this method is essentially 
identical to the Trotter-based adaptive time-dependent 
DMRG, with the advantages of finding the true varia- 
tional optimum and of global truncation after one entire 
infinitesimal time step as opposed to local truncations 
after each bond time step. Presently, it seems that the 
gain in precision is not extremely significant, and studies 
of relative computational costs for a desired precision are 
still lacking. 

6. Conclusion and Outlook 

I hope to have shown that there are now multiple ef- 
ficient ways of simulating Hamiltonian dynamics of pure 
and mixed quantum states for strongly correlated one- 
dimensional quantum systems in the DMRG framework. 



H+i- 



(51) 



This should open the way to a huge number of interesting 
new DMRG applications and be of interest for years to 
come. Very little is still known about the detailed perfor- 
mance of the new algorithms at the moment, but what 
is known is indeed very promising. In fact, I have passed 
over the fact that all the above can also extended to mas- 
ter equation dynamics of density matrices, i.e. dissipativc 
dynamics 8,9 - even less is known here at the moment. 
But algorithmically, everything is very close to the non- 
dissipative case: The approach of Zwolak and Vidal 9 is 
based on the observation that local (density) operators 
X) Pffff'l fJ )( fT 'l can be represented as N x N hermitian 
matrices. They now reinterpret the matrix coefficients as 
the coefficients of a local "superket" defined on a N 2 - 
dimensional local state space, and TEBD/DMRG-style 
dynamics can be carried out on that superket. Similarly, 
Verstraete et al. 8 consider dissipativc dynamics for states 
of the form of Eq. (48), which is formally equivalent to 
a matrix product state on a local N 2 dimensional state 
space. As should be clear from the discussion of matrix 
product states, both approaches are effectively identical 
but for the optimal approximation of Verstraete et al. 8 
replacing the somewhat less precise truncation of the ap- 
proach of Zwolak and Vidal. 9 Needless to say, even more 
interesting applications are waiting here! 
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